#
# More Donors, More Democracy
#
# Author: Sebastian Ziaja
# Date: 27 October 2018
#
# Replication script for
# Online Appendix C: Robustness checks
#
#

workdir <- ""
# workdir <- ifelse(basename(getwd()) == "aiddem_new", "", "../")
source(paste0(workdir, "scripts/load_packages.R"))
source(paste0(workdir, "scripts//extract.custom.R"))
source(paste0(workdir, "scripts//functions.R"))
source(paste0(workdir, "scripts//formulas.R"))
source(paste0(workdir, "scripts//variable_labels.R"))

load(paste0(workdir, "data/aiddem_datasets.RData"))
ddm <- read_rds(paste0(workdir, "data/temp/aiddyad_m.rds"))

iosh <- read_rds(paste0(workdir, "data/iosh.rds"))
iohq <- read_rds(paste0(workdir, "data/iohq.rds"))
mld <- read_csv(paste0(workdir, "data/multilateral_donors.csv"), col_types = "ccc")

dX  <- d1
dXm <- d1m

coefmap <- list(
            "l.CMgnh" = vln[2],
            "l.CMgap.log" = vla[2],
            "`l.CMgnh(fit)`" = vlne[2],
            "`I(l.CMgap.log * 10)(fit)`" = vlae[2],
            "`IA_l(fit)`" = vlies[2],
            "l.CMenh" = vln[3],
            "l.CMeap.log" = vla[3]
            # "l.pop.log.r" = vlctrls[1],
            # "l.gdpcap.log.r" = vlctrls[2],
            # "l.war25" = vlctrls[3]
        )


# Functions ---------------------------------------------

trg <- function(x, clevel = .95, ...)
{
        screenreg(x
                  , ci.force = TRUE
                  , ci.force.level = clevel
                  , digits = 2
                  , controls = TRUE
                  , ...)
}


# Base models -------------------------------

## Base specifications employed in this appendix --------------------------

runmodel(
    bf
    , data = dX %>% filter(smpl == 1)
) -> rb -> rbasecoef

trg(rb
    , custom.coef.map =
        list(
            "l.CMgnh" = vln[2],
            "l.CMgap.log" = vla[2],
            "`l.CMgnh(fit)`" = vlne[2],
            "`I(l.CMgap.log * 10)(fit)`" = vlae[2],
            "`IA_l(fit)`" = vlies[2],
            "l.CMenh" = vln[3],
            "l.CMeap.log" = vla[3],
            "l.pop.log.r" = vlctrls[1],
            "l.gdpcap.log.r" = vlctrls[2],
            "l.war25" = vlctrls[3]
        )
    )

## Base models excluding non-electoral regimes ------

runmodel(
    bf
    , data = dX %>% filter(smplR == 1)
) -> r -> rbasecoef

trg(r
    , custom.coef.map =
        list(
            "l.CMgnh" = vln[2],
            "l.CMgap.log" = vla[2],
            "`l.CMgnh(fit)`" = vlne[2],
            "`I(l.CMgap.log * 10)(fit)`" = vlae[2],
            "`IA_l(fit)`" = vlies[2],
            "l.CMenh" = vln[3],
            "l.CMeap.log" = vla[3],
            "l.pop.log.r" = vlctrls[1],
            "l.gdpcap.log.r" = vlctrls[2],
            "l.war25" = vlctrls[3]
        )
    )


## Full regression tables for models from the main text --------------

# OLS:

r.ols <- vector("list", 4)
r.ols[[1]] <- felm(v2x.polyarchy.n ~ l.CMgnh |
                    cnamef + periodf | 0 | cnamef
                , data = dX %>% filter(smpl == 1))
r.ols[[2]] <- felm(v2x.polyarchy.n ~ l.CMgnh +
                    l.CMgap.log + l.CMenh +
                    l.CMeap.log + l.pop.log.r + l.gdpcap.log.r + l.war25 |
                    cnamef + periodf | 0 | cnamef
                , data = dX %>% filter(smpl == 1))
r.ols[[3]] <- felm(v2x.polyarchy.n ~
                    l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25 |
                    cnamef + periodf | 0 | cnamef
                , data = dX %>% filter(smpl == 1))
r.ols[[4]] <- felm(v2x.polyarchy.n ~ IA_l + l.CMgnh +
                    l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25 |
                    cnamef + periodf | 0 | cnamef
                , data = dX %>% filter(smpl == 1))

screenreg(r.ols
, custom.coef.map =
    list(
            "l.CMgnh" = vln[2],
            "l.CMgap.log" = vla[2],
            "IA_l" = vlies[2],
            "l.CMenh" = vln[3],
            "l.CMeap.log" = vla[3],
            "l.pop.log.r" = vlctrls[1],
            "l.gdpcap.log.r" = vlctrls[2],
            "l.war25" = vlctrls[3]
    )
            , ci.force = TRUE
        , digits = 2
)

# First stage:

r1 <- vector("list", 4)

r1[[1]] <- felm(l.CMgnh ~ I(l.ZwvCMgwh94 / 10) |
                      cnamef + periodf | 0 |
                      cnamef, data = dX %>% filter(smpl == 1))
r1[[2]] <- felm(l.CMgnh ~ I(l.ZwvCMgwh94 / 10) +
                    l.pop.log.r + l.gdpcap.log.r + l.war25 |
                      cnamef + periodf | 0 |
                      cnamef, data = dX %>% filter(smpl == 1))
r1[[3]] <- felm(I(l.CMgap.log * 10) ~ I(l.ZwvCMgwh94 / 10) +
                    l.pop.log.r + l.gdpcap.log.r + l.war25
                    | cnamef + periodf | 0 |
                      cnamef, data = dX %>% filter(smpl == 1))
r1[[4]] <- felm(IA_l ~ I(l.ZwvCMgwh94 / 10) +
                    l.pop.log.r + l.gdpcap.log.r + l.war25
                    | cnamef + periodf | 0 |
                      cnamef, data = dX %>% filter(smpl == 1))

screenreg(r1
, custom.coef.map =
    list(
            "I(l.ZwvCMgwh94/10)" = vlzp[2],
            "l.pop.log.r" = vlctrls[1],
            "l.gdpcap.log.r" = vlctrls[2],
            "l.war25" = vlctrls[3]
    )
        , ci.force = TRUE
        , digits = 2
)

# Second stage:

r2 <- vector("list", 4)

r2[[1]] <- felm(v2x.polyarchy.n ~ 1 |
                      cnamef + periodf |
                      (l.CMgnh ~ l.ZwvCMgwh94) |
                      # (l.CMgnh ~ l.ZwvDistCM) |
                      cnamef, data = dX %>% filter(smpl == 1))
r2[[2]] <- felm(v2x.polyarchy.n ~
                    l.pop.log.r + l.gdpcap.log.r + l.war25 |
                      cnamef + periodf |
                      (l.CMgnh ~ l.ZwvCMgwh94) |
                      # (l.CMgnh ~ l.ZwvDistCM) |
                      cnamef, data = dX %>% filter(smpl == 1))
r2[[3]] <- felm(v2x.polyarchy.n ~ #1
                    # l.CMgnh +
                    l.pop.log.r + l.gdpcap.log.r + l.war25
                    | cnamef + periodf |
                      (I(l.CMgap.log * 1) ~ l.ZwvCMgwh94) |
                      cnamef, data = dX %>% filter(smpl == 1))
r2[[4]] <- felm(v2x.polyarchy.n ~ #1
                    # l.CMgnh +
                    l.pop.log.r + l.gdpcap.log.r + l.war25
                    | cnamef + periodf |
                      (IA_l ~ l.ZwvCMgwh94) |
                      cnamef, data = dX %>% filter(smpl == 1))

screenreg(r2,
       custom.coef.map = list(
            "`l.CMgnh(fit)`" = vlne[2],
            "`I(l.CMgap.log * 1)(fit)`" = vlae[2],
            "`IA_l(fit)`" = vlies[2],
            "l.pop.log.r" = vlctrls[1],
            "l.gdpcap.log.r" = vlctrls[2],
            "l.war25" = vlctrls[3]
       )
        , ci.force = TRUE
        , digits = 2
)


## Full regression tables for models from the main text excluding non-electoral regimes ------

# OLS

r.ols <- vector("list", 4)
r.ols[[1]] <- felm(v2x.polyarchy.n ~ l.CMgnh |
                    cnamef + periodf | 0 | cnamef
                , data = dX %>% filter(smplR == 1))
r.ols[[2]] <- felm(v2x.polyarchy.n ~ l.CMgnh +
                    l.CMgap.log + l.CMenh +
                    l.CMeap.log + l.pop.log.r + l.gdpcap.log.r + l.war25 |
                    cnamef + periodf | 0 | cnamef
                , data = dX %>% filter(smplR == 1))
r.ols[[3]] <- felm(v2x.polyarchy.n ~
                    l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25 |
                    cnamef + periodf | 0 | cnamef
                , data = dX %>% filter(smplR == 1))
r.ols[[4]] <- felm(v2x.polyarchy.n ~ IA_l + l.CMgnh +
                    l.CMgap.log + l.pop.log.r + l.gdpcap.log.r + l.war25 |
                    cnamef + periodf | 0 | cnamef
                , data = dX %>% filter(smplR == 1))

screenreg(r.ols
, custom.coef.map =
    list(
            "l.CMgnh" = vln[2],
            "l.CMgap.log" = vla[2],
            "IA_l" = vlies[2],
            "l.CMenh" = vln[3],
            "l.CMeap.log" = vla[3],
            "l.pop.log.r" = vlctrls[1],
            "l.gdpcap.log.r" = vlctrls[2],
            "l.war25" = vlctrls[3]
    )
            , ci.force = TRUE
        , digits = 2
)


# First stage

r1 <- vector("list", 4)

r1[[1]] <- felm(l.CMgnh ~ I(l.ZwvCMgwh94 / 10) |
                      cnamef + periodf | 0 |
                      cnamef, data = dX %>% filter(smplR == 1))
r1[[2]] <- felm(l.CMgnh ~ I(l.ZwvCMgwh94 / 10) +
                    l.pop.log.r + l.gdpcap.log.r + l.war25 |
                      cnamef + periodf | 0 |
                      cnamef, data = dX %>% filter(smplR == 1))
r1[[3]] <- felm(I(l.CMgap.log * 10) ~ I(l.ZwvCMgwh94 / 10) +
                    l.pop.log.r + l.gdpcap.log.r + l.war25
                    | cnamef + periodf | 0 |
                      cnamef, data = dX %>% filter(smplR == 1))
r1[[4]] <- felm(IA_l ~ I(l.ZwvCMgwh94 / 10) +
                    l.pop.log.r + l.gdpcap.log.r + l.war25
                    | cnamef + periodf | 0 |
                      cnamef, data = dX %>% filter(smplR == 1))

screenreg(r1
, custom.coef.map =
    list(
            "I(l.ZwvCMgwh94/10)" = vlzp[2],
            "l.pop.log.r" = vlctrls[1],
            "l.gdpcap.log.r" = vlctrls[2],
            "l.war25" = vlctrls[3]
    )
        , ci.force = TRUE
        , digits = 2
      )


# Second stage

r2 <- vector("list", 4)

r2[[1]] <- felm(v2x.polyarchy.n ~ 1 |
                      cnamef + periodf |
                      (l.CMgnh ~ l.ZwvCMgwh94) |
                      # (l.CMgnh ~ l.ZwvDistCM) |
                      cnamef, data = dX %>% filter(smplR == 1))
r2[[2]] <- felm(v2x.polyarchy.n ~
                    l.pop.log.r + l.gdpcap.log.r + l.war25 |
                      cnamef + periodf |
                      (l.CMgnh ~ l.ZwvCMgwh94) |
                      # (l.CMgnh ~ l.ZwvDistCM) |
                      cnamef, data = dX %>% filter(smplR == 1))
r2[[3]] <- felm(v2x.polyarchy.n ~ #1
                    # l.CMgnh +
                    l.pop.log.r + l.gdpcap.log.r + l.war25
                    | cnamef + periodf |
                      (I(l.CMgap.log * 1) ~ l.ZwvCMgwh94) |
                      cnamef, data = dX %>% filter(smplR == 1))
r2[[4]] <- felm(v2x.polyarchy.n ~ #1
                    # l.CMgnh +
                    l.pop.log.r + l.gdpcap.log.r + l.war25
                    | cnamef + periodf |
                      (IA_l ~ l.ZwvCMgwh94) |
                      cnamef, data = dX %>% filter(smplR == 1))

screenreg(r2,
       custom.coef.map = list(
            "`l.CMgnh(fit)`" = vlne[2],
            "`I(l.CMgap.log * 1)(fit)`" = vlae[2],
            "`IA_l(fit)`" = vlies[2],
            "l.pop.log.r" = vlctrls[1],
            "l.gdpcap.log.r" = vlctrls[2],
            "l.war25" = vlctrls[3]
       )
        , ci.force = TRUE
        , digits = 2
)


# Alternative indicators -----------------

## Aid flow types: commitments and disbursements --------------

xf <- bf %>%
    gsub("CM", "DM", .)
runmodel(xf, data = d3crs %>% filter(smpl == 1)) -> r

trg(r
, custom.coef.map = list(
            "l.DMgnh" = vln[2],
            "l.DMgap.log" = vla[2],
            "`l.DMgnh(fit)`" = vlne[2],
            "`I(l.DMgap.log * 10)(fit)`" = vlae[2],
            "`IA_l(fit)`" = vlies[2],
            "l.DMenh" = vln[3],
            "l.DMeap.log" = vla[3]
)
)

# excluding non-electoral regimes:

xf <- bf %>%
    gsub("CM", "DM", .)
runmodel(xf, data = d1crs %>% filter(smplR == 1)) -> r

trg(r
, custom.coef.map = list(
            "l.DMgnh" = vln[2],
            "l.DMgap.log" = vla[2],
            "`l.DMgnh(fit)`" = vlne[2],
            "`I(l.DMgap.log * 10)(fit)`" = vlae[2],
            "`IA_l(fit)`" = vlies[2],
            "l.DMenh" = vln[3],
            "l.DMeap.log" = vla[3]
)
    , omit.coef = "pop|gdp|war|X"
)

# both disbursement and commitment data

xf <-
list(paste0("v2x.polyarchy.n ~ l.CMgnh + l.DMgnh | cnamef + periodf | 0 | cnamef"),
     paste0("v2x.polyarchy.n ~ ",
            "l.CMgnh + l.CMenh + l.CMgap.log + l.CMeap.log +",
            "l.DMgnh + l.DMenh + l.DMgap.log + l.DMeap.log +",
            "l.pop.log.r + l.gdpcap.log.r + l.war25 | cnamef + periodf | 0 | cnamef")
     )

runmodel(xf, data = d1crs %>% filter(smpl == 1)) -> r

trg(r
    , custom.coef.names = c(paste0(vln[2], ", commitments"),
                            paste0(vln[2], ", disbursements"),
                            paste0(vln[3], ", commitments"),
                            paste0(vla[2], ", commitments"),
                            paste0(vla[3], ", commitments"),
                            paste0(vln[3], ", disbursements"),
                            paste0(vla[2], ", disbursements"),
                            paste0(vla[3], ", disbursements")
                        # , "pop", "gdp", "war"
                        )
    , reorder.coef = c(1, 2, 3, 6, 4, 7, 5, 8)
    , omit.coef = "pop|gdp|war"
)


## Aid donors: adding multilateral donors -------------------------------

runmodel(bf, data = dXm %>% filter(smpl == 1)) -> r

trg(r
    , custom.coef.map = coefmap
    , omit.coef = "pop|gdp|war|X"
)


# See "B1_spurious_placebo_test.R" for the script for the randomized placebo graphs


## Aid sectors: economic aid --------------------------

bx <-
        gsub("CMgnh |", "CMenh |",
        gsub("(l.CMgnh", "(l.CMenh",
        gsub("I(l.CMgap.log", "I(l.CMeap.log",
        gsub("l.ZwvCMgwh94", "l.ZlfCMewh94",
        gsub("IA_l", "IAe_l",
    bf
    , fixed = TRUE)
    , fixed = TRUE)
    , fixed = TRUE)
    , fixed = TRUE)
    , fixed = TRUE)


runmodel(bx
    , data = dX %>% filter(smpl == 1)) -> r

trg(r, custom.coef.map = list(
            "l.CMenh" = vln[3],
            "l.CMeap.log" = vla[3],
            "`l.CMenh(fit)`" = vlne[2],
            "`I(l.CMeap.log * 10)(fit)`" = vlae[2],
            "`IAe_l(fit)`" = vlies[2],
            "l.CMgnh" = vln[2],
            "l.CMgap.log" = vla[2]
        )
    )




## Aid sectors: relevant and irrelevant sub-sectors ---------------------------

# OLS models:

sectx <- c("g", "d", "c", "h", "t", "p")
xf <- vector("list", length(sectx))
bfx <- bf[[2]] %>%
        gsub("l.CMenh \\+ ", "", .) %>%
        gsub("l.CMeap.log \\+ ", "", .)
for(i in 1:length(sectx)) {
    bfx %>% gsub("CMg", paste0("CM", sectx[i]), .) -> xf[[i]]
}

runmodel(xf) -> rb

trg(
    rb
    , custom.model.names = c("Democracy", "Democracy, no PB", "Civil society",
                             "Health", "Transport", "Energy")
    , custom.coef.names = c(rep("No. of donors", 6))
    , omit.coef = "log|war"
)

# IV models:

sectx <- c("g", "d", "c", "h", "t", "p")
xf <- vector("list", length(sectx))
bfx <- bf[[3]] %>%
        gsub("l.CMenh \\+ ", "", .) %>%
        gsub("l.CMeap.log \\+ ", "", .)
for(i in 1:length(sectx)) {
    bfx %>% gsub("CMg", paste0("CM", sectx[i]), .) -> xf[[i]]
}

runmodel(xf, data = dX %>% filter(smpl == 1)) -> rb

trg(
    rb
    , custom.model.names = c("Democracy", "Democracy, no PB", "Civil society",
                             "Health", "Transport", "Energy")
    , custom.coef.names = c(rep("No. of donors", 6))
    , omit.coef = "log|war"
)

# IV models excluding non-electoral systems:

runmodel(xf, data = dX %>% filter(smplR == 1)) -> rb

trg(rb
    , custom.model.names = c("Democracy", "Democracy, no PB", "Civil society",
                             "Health", "Transport", "Energy")
    , custom.coef.names = c(rep("No. of donors", 6))
    , omit.coef = "log|war"
)


# Effect on liberal and participatory components:

sectx <- c("g", "d", "c", "g", "d", "c")
xf <- vector("list", length(sectx))
bfx <- bf[[3]] %>%
        gsub("l.CMenh \\+ ", "", .) %>%
        gsub("l.CMeap.log \\+ ", "", .)
for(i in 1:length(sectx)) {
    bfx %>% gsub("CMg", paste0("CM", sectx[i]), .) -> xf[[i]]
}
for(i in 1:3) {
    xf[[i]] %>% gsub("v2x.polyarchy.n", "v2x_liberal_n", .) -> xf[[i]]
}
for(i in 4:6) {
    xf[[i]] %>% gsub("v2x.polyarchy.n", "v2x_partip_n", .) -> xf[[i]]
}

runmodel(xf, data = dX %>% filter(smpl == 1)) -> rb

trg(
    rb
    , custom.coef.names = c(rep(c("Democracy", "Democracy, no PB", "Civil society"), 1))
    , omit.coef = "log|war"
)

## Aid amount references: per GDP and per capita ------------------------

bf %>%
    gsub("gap.log", "ggp.log", .) %>%
    gsub("eap.log", "egp.log", .) -> xf
bf %>%
    gsub("gap.log", "gcp.log", .) %>%
    gsub("eap.log", "ecp.log", .) -> xf2
xf[[2]] -> xf[[1]]
xf2[[2]] -> xf[[2]]
xf[[4]] -> xf[[3]]
xf2[[4]] -> xf[[4]]
sub("IA", "IApc", bf[[5]]) -> xf[[5]]
sub("IA", "IAgdp", bf[[5]]) -> xf[[6]]

runmodel(xf) -> r

screenreg(r
       , custom.coef.map = list(
           "l.CMgnh" = vln[2],
           "`l.CMgnh(fit)`" = vlne[2],
           "l.CMggp.log" = sub("USD per capita", "share of GDP", vlc[2]),
           "`I(l.CMggp.log * 10)(fit)`" = "Dem. aid (share of GDP) * 10",
           "l.CMgcp.log" = vlc[2],
           "`I(l.CMgcp.log * 10)(fit)`" = paste(vlc[2], "* 10"),
           "`IApc_l(fit)`" = "No. dem. don.  *  Dem. aid p.c. (est.)",
           "`IAgdp_l(fit)`" = "No. dem. don.  *  Dem. aid GDP (est.)",
           "l.CMenh" = vln[3],
           "l.CMegp.log" = sub("USD per capita", "share GDP", vlc[3]),
           "l.CMecp.log" = vlc[3]

)
           , omit.coef = "pop|gdp|war"
       , digits = 2
    , ci.force = TRUE
)



## Aid fragmentation indicators ------------------------------

# Small donors:

bf %>%
    sub("l.CMgnh", "l.CMgsp", .) %>%
    sub("l.CMenh", "l.CMesp", .) %>%
    sub("IA_l", "IAsmall", .) -> xf

runmodel(xf, data = d1frag %>% filter(smpl == 1)) -> r

trg(r
    , custom.coef.map = list(
            "l.CMgsp" = vln[2] %>% sub("of", "of small", .),
            "l.CMgap.log" = vla[2],
            "`l.CMgsp(fit)`" = vlne[2] %>% sub("of", "of small", .),
            "`I(l.CMgap.log * 10)(fit)`" = vlae[2],
            "`IAsmall(fit)`" = vlies[2] %>% sub("of", "of small", .),
            "l.CMesp" = vln[3] %>% sub("of", "of small", .),
            "l.CMeap.log" = vla[3]
        )
, omit.coef = "pop|gdp|war|X"
)

# Excluding non-electoral regimes:

bf %>%
    sub("l.CMgnh", "l.CMgsp", .) %>%
    sub("l.CMenh", "l.CMesp", .) %>%
    sub("IA_l", "IAsmall", .) -> xf

runmodel(xf, data = d1frag %>% filter(smplR == 1)) -> r

trg(r
    , custom.coef.map = list(
            "l.CMgsp" = vln[2] %>% sub("of", "of small", .),
            "l.CMgap.log" = vla[2],
            "`l.CMgsp(fit)`" = vlne[2] %>% sub("of", "of small", .),
            "`I(l.CMgap.log * 10)(fit)`" = vlae[2],
            "`IAsmall(fit)`" = vlies[2] %>% sub("of", "of small", .),
            "l.CMesp" = vln[3] %>% sub("of", "of small", .),
            "l.CMeap.log" = vla[3]
        )
    , omit.coef = "pop|gdp|war|X"
)



# inverted HHI

bf %>%
    gsub("l\\.CMgnh", "l.CMgfp", .) %>%
    gsub("l\\.CMefp", "l.CMenp", .) %>%
    sub("IA_l", "IAhhi", .) -> xf
xf[[3]] %>% gsub("l\\.CMgfp", "I(l.CMgfp * 100)", .) -> xf[[3]]
xf[[4]] %>% gsub("l\\.CMgfp", "I(l.CMgfp * 100)", .) -> xf[[4]]

runmodel(xf, data = d1frag %>% filter(smplR == 1)) -> r

trg(r
, custom.coef.map =
    list(
            "l.CMgfp" = vlf[2],
            "l.CMgap.log" = vla[2],
            "`I(l.CMgfp * 100)(fit)`" = vlf[2] %>% paste(., "(est.)"),
            "`I(l.CMgap.log * 10)(fit)`" = vlae[2],
            "`IAhhi(fit)`" = paste(vlfs[2], " * ", vlas[2], "(est.)"),
            "l.CMefp" = vlf[3],
            "l.CMeap.log" = vla[3]
            # "l.pop.log.r" = vlctrls[1],
            # "l.gdpcap.log.r" = vlctrls[2],
            # "l.war25" = vlctrls[3]
        )
    , omit.coef = "pop|gdp|war|X"
)

# reversed concentration ratio.

bf %>%
    sub("l.CMgnh", "l.CMgrp", .) %>%
    sub("l.CMefp", "l.CMerp", .) %>%
    sub("IA_l", "IAcr3", .) -> xf
xf[[3]] %>% gsub("l\\.CMgrp", "I(l.CMgrp * 100)", .) -> xf[[3]]
xf[[4]] %>% gsub("l\\.CMgrp", "I(l.CMgrp * 100)", .) -> xf[[4]]

runmodel(xf, data = d1frag %>% filter(smpl == 1)) -> r

trg(r
    , custom.coef.map =
    list(
            "l.CMgrp" = vlf[2] %>% gsub("fragmentation", "concentration (rev.)", .),
            "l.CMgap.log" = vla[2],
            "`I(l.CMgrp * 100)(fit)`" = vlf[2] %>% paste(., "(est.)") %>% gsub("fragmentation", "concentration (rev.)", .),
            "`I(l.CMgap.log * 10)(fit)`" = vlae[2],
            "`IAcr3(fit)`" = paste(vlfs[2], " * ", vlas[2], "(est.)") %>% gsub("frag", "concentr. rev", .),
            "l.CMefp" = vlf[3] %>% gsub("fragmentation", "concentration (rev.)", .),
            "l.CMeap.log" = vla[3]
            # "l.pop.log.r" = vlctrls[1],
            # "l.gdpcap.log.r" = vlctrls[2],
            # "l.war25" = vlctrls[3]
        )
)


## Democracy indicators ----------------------------

# Polity:

bf %>%
    gsub("v2x.polyarchy.n", "e_polity2", .) -> xf

runmodel(xf) -> r

trg(r
    , custom.coef.map = coefmap
    , omit.coef = "pop|gdp|war|X"
)

# Freedom House:

bf %>%
    gsub("v2x.polyarchy.n", "fh.i.n", .) -> xf

runmodel(xf) -> r

trg(r
    , custom.coef.map = coefmap
    , omit.coef = "pop|gdp|war|X"
)

# Unified Democracy Scale:

bf %>%
    gsub("v2x.polyarchy.n", "uds.n", .) -> xf

runmodel(xf) -> r

trg(r
    , custom.coef.map = coefmap
    , omit.coef = "pop|gdp|war|X"
)


## Sub-components of democracy -----------------

xf <- vector("list", 1)
subx <- c("v2x_liberal_n", "v2x_partip_n", "v2xdl_delib_n", "v2x_egal_n")
for (i in 1:length(subx)) {
    xf[[i]] <- sub("v2x.polyarchy.n", subx[i], bf[[2]])
}

runmodel(xf) -> r

trg(r
    , custom.coef.names = vlsets[-(5:8)]
    , omit.coef = "pop|gdp|war"
    , custom.model.names = c("Liberal", "Participatory", "Deliberative", "Egalitarian")
)

xf <- vector("list", 1)
for (i in 1:length(subx)) {
    xf[[i]] <- sub("v2x.polyarchy.n", subx[i], bf[[3]])
}

runmodel(xf) -> r

trg(r
    , custom.coef.map = list("`l.CMgnh(fit)`" = vlsets[8])
    , custom.model.names = c("Liberal", "Participatory", "Deliberative", "Egalitarian")
)


## Indicators of democratic diversity ---------------------------

xf <- vector("list", 1)
subx <- c("v2xsharelargestparty", "v2xnopartyplatforms",
          "v2xrespectarguments", "v2xengagedsociety",
          "v2xcivsocpartenviron")
for (i in 1:length(subx)) {
    xf[[i]] <- sub("v2x.polyarchy.n", paste0("I(", subx[i], " * 10)"), bf[[2]])
}

runmodel(xf) -> r

trg(r
    , custom.coef.names = vlsets[-(5:8)]
    , omit.coef = "pop|gdp|war"
    , custom.model.names = c("Largest party share", "Party platforms",
                             "Counterarguments", "Engaged society",
                             "CSO environment")
)

xf <- vector("list", 1)
for (i in 1:length(subx)) {
    xf[[i]] <- sub("v2x.polyarchy.n", paste0("I(", subx[i], " * 10)"), bf[[3]])
}

runmodel(xf) -> r

trg(r
    , custom.coef.map = list("`l.CMgnh(fit)`" = vlsets[8])
    , custom.model.names = c("Largest party share", "Party platforms",
                             "Counterarguments", "Engaged society",
                             "CSO environment")
)

## Additional controls ----------------------------------


vector("list", 3) -> xf

xf <- NULL
bf[[2]] %>% sub("l.war25",
                "l.war25 + I(l.tradegdp / 10) + I(l.gdpgrowth / 1) + I(l.eu.prospect / 1)", .) -> xf[[1]]
xf[[1]] %>% sub("prospect / 1)",
                "prospect / 1) + I(l.fuelexports / 1)", .) -> xf[[1]]
bf[[3]] %>%
    gsub("1", "l.pop.log.r + l.gdpcap.log.r + l.war25", .) %>%
    sub("l.war25",
                "l.war25 + I(l.tradegdp / 10) + I(l.gdpgrowth / 1) + I(l.eu.prospect / 1)", .) -> xf[[2]]
xf[[2]] %>% sub("prospect / 1)",
                "prospect / 1) + I(l.fuelexports / 1)", .) -> xf[[2]]
xf[[2]] %>% sub("l.CMgnh", "IA_l", .) -> xf[[3]]

runmodel(xf) -> r

trg(r
    , custom.coef.map =
        list(
            "l.CMgnh" = vln[2],
            "l.CMgap.log" = vla[2],
            "`l.CMgnh(fit)`" = vlne[2],
            # "`I(l.CMgap.log * 10)(fit)`" = vlae[2],
            "`IA_l(fit)`" = vlies[2],
            "l.CMenh" = vln[3],
            "l.CMeap.log" = vla[3],
            "l.pop.log.r" = vlctrls[1],
            "l.gdpcap.log.r" = vlctrls[2],
            "l.war25" = vlctrls[3],
            "I(l.tradegdp/10)" = "Trade (percent of GDP) x 1/10",
                "I(l.gdpgrowth/1)" = "GDP growth (percent)",
            "I(l.eu.prospect/1)" = "EU accession prospect",
            "I(l.fuelexports/1)" = "Fuel exports (percent of exports)"
        )
)


# Alternative samples -------------------
## Alternative temporal units ------------------------

# Two-year periods:

bf -> xf
runmodel(xf, data = d2 %>% filter(smpl == 1)) -> r

trg(r
, custom.coef.map = coefmap
)


# Three-year periods:

bf -> xf
runmodel(xf, data = d3 %>% filter(smpl == 1)) -> r

trg(r
, custom.coef.map = coefmap
)

# excluding non-electoral regimes:

runmodel(xf, data = d3 %>% filter(smplR == 1)) -> r

trg(r
, custom.coef.map = coefmap
)

# Four-year periods:

bf -> xf
runmodel(xf, data = d4 %>% filter(smpl == 1)) -> r

trg(r
, custom.coef.map = coefmap
    , omit.coef = "pop|gdp|war|X"
)


## Temporal sub-samples -----------------------------------------

# 1994-2003

runmodel(bf, data = d1 %>% filter(smpl == 1, year < 2004)) -> r

trg(r
, custom.coef.map = coefmap
    , omit.coef = "pop|gdp|war|X"
)

# 2004-2013

runmodel(bf, data = d1 %>% filter(smpl == 1, year > 2003
                                  )) -> r

trg(r
, custom.coef.map = coefmap
    , omit.coef = "pop|gdp|war|X"
)

## Geographic sub-samples -----------------------------

# SSA:

bf -> xf
runmodel(xf, data = dX %>% filter(smpl == 1,
    region_ssa == 1)) -> r

trg(r
, custom.coef.map = coefmap
    , omit.coef = "pop|gdp|war|X"
)

# Not SSA:

bf -> xf
runmodel(xf, data = dX %>% filter(smpl == 1, region_ssa == 0)) -> r

trg(r
, custom.coef.map = coefmap
    , omit.coef = "pop|gdp|war|X"
)

# Not SSA, excluding non-electoral regimes:

bf -> xf
runmodel(xf, data = dX %>% filter(smplR == 1, region_ssa == 0)) -> r

trg(r
, custom.coef.map = coefmap
    , omit.coef = "pop|gdp|war|X"
)


# Alternative models ---------------------------------

## Lagged dependent variable (LDV) models ----------------------------

# Without controls:

xf <-
    bf[[1]] %>%
    gsub("cnamef \\+ periodf", "periodf", .) %>%
    gsub("~", "~ l.v2x.polyarchy.n +", .)
r <- NULL
felm(as.formula(xf), data = d1 %>% filter(smpl == 1)) -> r[[1]]
felm(as.formula(xf), data = d2 %>% filter(smpl == 1)) -> r[[2]]
felm(as.formula(xf), data = d3 %>% filter(smpl == 1)) -> r[[3]]

trg(r , custom.coef.map = list(
     "(Intercept)" = "(Intercept)",
     "l.v2x.polyarchy.n" = "Electoral democracy",
            "l.CMgnh" = vln[2])
    , omit.coef = "pop|gdp|war|X"
    , custom.model.names = c("One-year periods", "Two-year periods", "Three-year periods")
)

# With controls:

xf <-
    bf[[2]] %>%
    gsub("cnamef \\+ periodf", "periodf", .) %>%
    gsub("~", "~ l.v2x.polyarchy.n +", .)
r <- NULL
felm(as.formula(xf), data = d1 %>% filter(smpl == 1)) -> r[[1]]
felm(as.formula(xf), data = d2 %>% filter(smpl == 1)) -> r[[2]]
felm(as.formula(xf), data = d3 %>% filter(smpl == 1)) -> r[[3]]

trg(r
 , custom.coef.map = list(
     "(Intercept)" = "(Intercept)",
     "l.v2x.polyarchy.n" = "Electoral democracy",
            "l.CMgnh" = vln[2],
            "l.CMenh" = vln[3],
            "l.CMgap.log" = vla[2],
            "l.CMeap.log" = vla[3]
        )
    , custom.model.names = c("One-year periods", "Two-year periods", "Three-year periods")
)


## Nonlinear effects ---------------------

vector("list", 2) -> xf -> r

bf[[1]] -> xf[[1]]
bf[[2]] -> xf[[2]]
xf %>% gsub("l.CMgnh",
                 "l.CMgnh + I((l.CMgnh^2) / 100)", .) -> xf
runmodel(xf) -> r

trg(r
    , custom.coef.names = c(vlns[2], paste0(vlns[2], " (squared) x 1/100"),
                            vlns[3], vlas[2], vlas[3])
    , omit.coef = "pop|gdp|war"
)

# Interaction plot:

paste0("v2x.polyarchy.n ~ l.CMgnh + I(l.CMgnh^2)",
       " + l.CMenh + l.CMgap.log + l.CMeap.log + l.pop.log.r",
       " + l.gdpcap.log.r + l.war25",
       " + cnamef + periodf") -> xfp
lm(xfp, data = dX %>% filter(smpl == 1)) -> rp

interplot(rp, var1 = "l.CMgnh", var2 = "l.CMgnh",
          sims = 5000, plot = FALSE, hist = TRUE) -> ip
# get ratio of non-clustered to clustered coefs:
summary(r[[2]])$coefficients[1, 2] / (1 * summary(rp)$coefficients[2, 2]) -> iprate
# correct se by ratio:
ip$ub <- ip$coef + ((ip$ub - ip$coef) * iprate)
ip$lb <- ip$coef - ((ip$coef - ip$lb) * iprate)

names(ip)[1:2] <- c("fake", "coef1")  # bug in interplot(); must rename

interplot(ip, hist = TRUE, var2_dt = dX %>% filter(smpl == 1) %>% select(l.CMgnh) %>% unlist) +
    xlab('Number of democracy donors') +
    ylab('Estimated coefficient for number of democracy donors') +
    theme_bw() +
    geom_hline(yintercept = 0, linetype = 2)


## Interactions with democracy ------------------------------------

r1 <- vector("list", 1)

r1[[1]] <- felm(v2x.polyarchy.n ~ l.CMgnh + l.v2x.polyarchy.n  +
                    I(l.CMgnh  * l.v2x.polyarchy.n  / 10)
       | periodf
       | 0 | cnamef, data = dX %>% filter(smpl == 1))
r1[[2]] <- felm(v2x.polyarchy.n ~ l.CMgnh + l.v2x.polyarchy.n +
                    I(l.CMgnh  * l.v2x.polyarchy.n  / 10) +
                    l.pop.log.r + l.gdpcap.log.r + l.war25 | periodf
                | 0 |
                    cnamef, data = dX %>% filter(smpl == 1))
r1[[3]] <- felm(v2x.polyarchy.n ~ l.CMgnh + l.v2x.polyarchy.n +
                    I(l.CMgnh  * l.v2x.polyarchy.n  / 10) + I(l.CMgap.log * 1) +
                    l.pop.log.r + l.gdpcap.log.r + l.war25 | periodf
                | 0 |
                    cnamef, data = dX %>% filter(smpl == 1))


trg(r1
    , custom.coef.map = list(
        "l.v2x.polyarchy.n" = "Electoral democracy",
            "l.CMgnh" = vln[2],
        "I(l.CMgnh * l.v2x.polyarchy.n/10)" = paste(vlns[2], "* Elect. dem."),
            "I(l.CMgap.log * 1)" = vla[2]
            # "l.pop.log.r" = vlctrls[1],
            # "l.gdpcap.log.r" = vlctrls[2],
            # "l.war25" = vlctrls[3]
        )
    , omit.coef = "period|cname"
)

# Interaction plot:

paste0("v2x.polyarchy.n ~ 1 + l.CMgnh * l.v2x.polyarchy.n",
       " + l.pop.log.r + l.gdpcap.log.r + l.war25"
       , " + periodf"
       ) -> xfp
lm(xfp, data = dX %>% filter(smpl == 1)) -> rp

rtemp <- felm(v2x.polyarchy.n ~
                    l.CMgnh  * l.v2x.polyarchy.n + #l.CMgap.log +
                    l.pop.log.r + l.gdpcap.log.r + l.war25 | periodf
                | 0 |
                    cnamef, data = dX %>% filter(smpl == 1))

# cluster se at country level
interplot(rp, var1 = "l.CMgnh", var2 = "l.v2x.polyarchy.n",
          sims = 5000, plot = FALSE, hist = TRUE) -> ip
# get ratio of non-clustered to clustered coefs:
summary(rtemp)$coefficients[1, 2] / (1 * summary(rp)$coefficients[2, 2]) -> iprate
# correct se by ratio:
ip$ub <- ip$coef + ((ip$ub - ip$coef) * iprate)
ip$lb <- ip$coef - ((ip$coef - ip$lb) * iprate)

names(ip)[1:2] <- c("fake", "coef1")  # bug in interplot(); must rename

interplot(ip, hist = TRUE, var2_dt = dX$l.v2x.polyarchy.n[dX$smpl == 1]) +
    xlab('Democracy (V-Dem polyarchy, rescaled to [0; 100]') +
    ylab('Estimated coefficient for number of democracy donors') +
    theme_bw() +
    geom_hline(yintercept = 0, linetype = 2)


## Time trends -----------------------------------------

# Global:

bf %>%
    gsub("~", "~ period +", .) %>% gsub("cnamef \\+ periodf", "cnamef", .) -> xf

runmodel(xf, kcx = TRUE) -> r

trg(r
    , custom.coef.map = c(list("period" = "Time trend"), coefmap)
)

# Regional:

dX %>% filter(reg == "..") %>% select(cc) %>% arrange(cc) %>%
    unique %>% c %>% unlist %>% sort -> clist

for(i in c(3, 4, 6, 7, 9, 11)) {
    dX$reg[dX$cc == clist[i]] <- "EuropeCentralAsia"
}

for(i in c(8, 10)) {
    dX$reg[dX$cc == clist[i]] <- "MiddleEastNorthAfrica"
}

dX$reg[dX$cc == clist[1]] <- "EastAsiaPacific"
dX$reg[dX$cc == clist[2]] <- "LatinAmericaCaribbean"
dX$reg[dX$cc == clist[5]] <- "SubSaharanAfrica"
dX$reg[dX$cc == clist[12]] <- "LatinAmericaCaribbean"

regl <- unique(dX$reg)[-2]
for(i in 1:length(regl)) {
    dX[, paste0("trend", regl[i])] <- dX$period * as.numeric(dX$reg == regl[i])
}

bf %>%
    gsub("~",
         paste0("~ ",
                paste(grep("trend", names(dX), value = TRUE), collapse = " + "),
                " +"), .) -> xf

runmodel(xf) -> r

trg(r
    , custom.coef.map =
        c(
            list("trendSouthAsia" = "Time trend, South Asia"),
            list("trendEuropeCentralAsia" = "Time trend, Europe and Central Asia"),
            list("trendLatinAmericaCaribbean" = "Time trend, Latin America and Caribbean"),
            list("trendEastAsiaPacific" = "Time trend, East Asia and Pacific"),
            list("trendMiddleEastNorthAfrica" = "Time trend, Middle East and North Africa"),
            coefmap
        )
    , omit.coef = "pop|gdp|war"
)


## Contemporaneous effects -----------------------------

bf %>%
    gsub("l\\.", "", .) -> xf
runmodel(xf, data = dX %>% filter(smpl == 1)) -> r

coefmap2 <- coefmap
names(coefmap2) <- names(coefmap) %>% gsub("l.", "", ., fixed = TRUE)

trg(r
    , custom.coef.map = coefmap2
    , omit.coef = "pop|gdp|war|X"
)

# excluding non-electoral regimes:

bf %>%
    gsub("l\\.", "", .) -> xf
runmodel(xf, data = dX %>% filter(smplR == 1)) -> r

coefmap2 <- coefmap
names(coefmap2) <- names(coefmap) %>% gsub("l.", "", ., fixed = TRUE)

trg(r
    , custom.coef.map = coefmap2
    , omit.coef = "pop|gdp|war|X"
)


## Two-way clustered standard errors -----------------------------

bbf <- vector("list", 1)

bbf0 <- "v2x.polyarchy.n ~ 1 | cnamef + periodf | 0 | cnamef + periodf"
bbf[[1]] <- sub("1", "l.CMgnh", bbf0)
bbf[[2]] <- sub("1",
               paste0("l.CMgnh + l.CMenh + l.CMgap.log + ",
                      "l.CMeap.log + l.pop.log.r + l.gdpcap.log.r + l.war25"),
                      bbf0)
bbf[[3]] <- sub("0",
               "(l.CMgnh ~ l.ZwvCMgwh94)",
               bbf0)
bbf[[4]] <- sub("0",
               "(I(l.CMgap.log * 10) ~ l.ZwvCMgwh94)",
                      bbf0)
bbf[[5]] <- sub("0",
               "(IA_l ~ l.ZwvCMgwh94)",
                      bbf0)

runmodel(bbf, data = d1 %>% filter(smpl == 1)) -> r

trg(r
    , custom.coef.map = coefmap
    , omit.coef = "pop|gdp|war"
    , caption = "Base models, clustered by country and time period"
)


## Heteroscedasticity-and-autocorrelation-robust standard errors -------------

# See "C1_conley_se.R" for calculation

load("data//conleyin_d1.RData")

screenreg(r
    , custom.coef.map = coefmap
       , omit.coef = "pop|gdp|war"
       , override.ci.low = haccilow
       , override.ci.up = haccihi
       , digits = 2
)


## Reverse models ----------------------------

r.rev <- vector("list", 1)
r.rev[[1]] <- felm(CMgnh ~ I(l.v2x.polyarchy.n / 10) |
                    cnamef + periodf | 0 | cnamef
                , data = dX %>% filter(smpl == 1))
r.rev[[2]] <- felm(CMgnh ~ I(l.v2x.polyarchy.n / 10) +
                    l.pop.log.r + l.gdpcap.log.r + l.war25 |
                    cnamef + periodf | 0 | cnamef
                , data = dX %>% filter(smpl == 1))
r.rev[[3]] <- felm(CMgnh ~ I(l.v2x.polyarchy.n / 10) + l.CMenh +
                    l.CMgap.log + l.CMeap.log + l.pop.log.r + l.gdpcap.log.r + l.war25 |
                    cnamef + periodf | 0 | cnamef
                , data = dX %>% filter(smpl == 1))

screenreg(r.rev
    , custom.coef.names = c("Level of democracy / 10", vla[2], vlctrls, vln[3], vla[3])
    , reorder.coef = c(1:2, 6:7, 3:5)
        , ci.force = TRUE
        , digits = 2
)

